//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   SweepSource.hh
 * \author Jeremy Roberts
 * \date   Aug 18, 2011
 * \brief  Right hand side for sweeping.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev:: 126                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-04 23:13:48 +0000 (Sun, 04 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#ifndef SWEEPSOURCE_HH_
#define SWEEPSOURCE_HH_

#include "angle/AngularMesh.hh"
#include "angle/MomentsIndex.hh"
#include "angle/Moment_to_Discrete.hh"
#include "angle/Octant.hh"

#include "tran/Definitions.hh"
#include "tran/State.hh"
#include "tran/CartMesh.hh"
#include "tran/Moments_Field.hh"
#include "tran/Psi_Cell_Field.hh"

#include "material/Materials.hh"

#include "source/BaseSource.hh"
#include "tran/FissionSource.hh"

#include "../utilities/SP.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class SweepSource
 * \brief Class for sweep sources.
 *
 * This class defines a general right-hand-side construct for sweeps (which
 * are the basis for all solvers).  That is, the SweepSource is a source
 * for a sweep along a particular angle.
 *
 * The source for a given sweep is comprised in general of these four
 * components:
 *  - within-group scattering source
 *  - in-scattering source (from both up- and down-scattering)
 *  - fission source
 *  - external source
 * Algorithmically, the fission source behaves as an external source within
 * an inner iteration.  Moreover, and of importance for response function
 * generation, there \em can be a fission source (i.e. multiplication)
 * together with an fixed source.
 *
 * Recall that a within group equation is represented as
 * \f[
 *    \mathbf{T}[\psi]_g = [\mathbf{M}][\mathbf{S}]_{gg}[\phi]_g + \bar{Q}_g \, ,
 * \f]
 * where \f$ \bar{Q}_g \f$ represents everything \em but the within-group
 * scattering.  The right hand side is the sweep source (but \em not the
 * right hand side for the linear system of interest, which casts the
 * problem in terms of flux moments!).
 *
 * As a sanity check, consider a purely
 * isotropic flux in one group with isotropic scattering. Then we must
 * have the equation
 * \f[
 *    \mathbf{T}\psi = \Sigma_{0} \phi_{00} / 4\pi \, ,
 * \f]
 * which indicates that application of \f$ \mathbf{M} \f$ implies
 * a normalization.
 *
 * For each inner iteration, the scattering components are stored in
 * moments form.  The M operator is then applied using one row at a
 * time for the associated angle being swept.
 *
 * \note Keep in mind where we're headed: solving the linear system
 *       \f$ (\mathbf{I}-\mathbf{D}\mathbf{T}^{-1}\mathbf{M}\mathbf{S})\phi
 *        = \mathbf{D}\mathbf{T}^{-1}q \f$. Hence, the sweep source is q,
 *        and we'll invert the transport operator \em T via sweeps.
 *
 * \note Denovo uses a different class for each source component, and then
 *       a "database" to house all the relevant components.  Here, we'll
 *       keep everything in a single class, with different methods that
 *       will set various components.  Of course, fission and external
 *       sources are still defined externally.
 *
 * \sa InScatterSource, WithinGroupSource, FissionSource, BaseSource
 */
/*!
 * \example test/testSweepSource.cc
 *
 * Test of SweepSource.
 */
//===========================================================================//
template <class Equations, class D>
class SweepSource
{

public:

  /// Typedefs
  //\{
  //
  typedef typename util::SP<SweepSource<Equations, D> >             SP_sweepsource;
  //
  typedef typename State<Equations, D>::SP_state                    SP_state;
  typedef typename Base_Mesh::SP_mesh                               SP_mesh;
  typedef typename AngularMesh<D>::SP_angularmesh                   SP_angularmesh;
  typedef typename Octant<D>::Omega                                 Omega;
  typedef Materials::SP_materials                                   SP_materials;
  //
  typedef typename BaseSource<D>::SP_source                         SP_externalsource;
  typedef typename FissionSource<Equations, D>::SP_fissionsource    SP_fissionsource;
  //
  typedef Moments_Field<Equations, D>                               Moments_Field_t;
  typedef Psi_Cell_Field<Equations, D>                              Discrete_Field_t;
  typedef Cell_Field<double>                                        SweepSourceField;
  //
  typedef typename Moment_to_Discrete<D>::SP_moment_to_discrete     SP_moment_to_discrete;
  typedef typename Moment_to_Discrete<D>::M_Row                     M_Row;
  typedef typename MomentsIndex<D>::SP_momentsindex                 SP_momentsindex;
  //\}

  /*!
   * \brief Constructor.
   *
   * This sizes the source variables.
   *
   * \param state           The state.
   * \param mesh            The mesh.
   * \param angularmesh     The angular mesh.
   * \param materials       The material library.
   * \param momentsindex    The moments index.
   * \param m_operator      The moments to discrete operator.
   *
   */
  SweepSource(SP_state state,
              SP_mesh mesh,
              SP_angularmesh angularmesh,
              SP_materials materials,
              SP_momentsindex momentsindex,
              SP_moment_to_discrete m_operator);

  /*!
   * \brief Set an external source.
   *
   * \param source  Smart pointer to external source
   */
  void set_external_source(SP_externalsource source)
  {
    d_externalsource = source;
  }

  /*!
   * \brief Set a fission source.
   *
   * \param source  Smart pointer to fission source
   */
  void set_fission_source(SP_fissionsource source)
  {
    d_fissionsource = source;
  }

  /*!
   * \brief Build the in-scatter source in terms of moments.
   *
   * \param g       Group of problem we are solving.
   * \param phi_g   Iteration flux moments vector.
   */
  void build_in_scatter_source(int g);

  /*!
   * \brief Build the within-group source in terms of moments.
   *
   * \param g       Group of problem we are solving.
   * \param phi_g   Iteration flux moments vector.
   */
  void build_within_group_source(int g, const Moments_Field_t &phi_g);

  /*!
   * \brief Build the complete sweep source for this group.
   *
   * \param g       Group of problem we are solving.
   * \param o       Octant of angle we are sweeping.
   * \param a       Angle we are sweeping.
   *
   */
   void build_source(int g, int o, int a);

   /// \name Accessors
   //\{

   /// Return sweep source.
   const SweepSourceField& source() const
   {
     return d_source;
   }

   /// Return sweep source.
   const Moments_Field_t& within_group_source() const
   {
     return d_within_group_source;
   }

   /// Return sweep source.
   const Moments_Field_t& in_scatter_source() const
   {
     return d_in_scatter_source;
   }

   /// Return state.
   SP_state state() const
   {
     return d_state;
   }

   /// Return mesh.
   SP_mesh mesh() const
   {
     return d_mesh;
   }

   /// Return moments index.
   SP_momentsindex momentsindex() const
   {
     return d_momentsindex;
   }

   //\}

private:

  // >>> Data

  //---------------------------------------------------------------------------
  // Local source data

  /// Total sweep source over all unknowns for a single angle.
  SweepSourceField d_source;

  /// Within-group scattering source over all unknowns in moment form.
  Moments_Field_t d_within_group_source;

  /// In-scatter scattering source over all unknowns in moment form.
  Moments_Field_t d_in_scatter_source;

  //---------------------------------------------------------------------------
  // Source group and angle.

  /// Group.
  int d_group;

  /// Octant index.
  int d_octant;

  /// Angle index in octant.
  int d_angle;


  //---------------------------------------------------------------------------
  // Smart pointers to external/fission sources.

  /// Pointer to external source.
  SP_externalsource d_externalsource;

  /// Pointer to fission source.
  SP_fissionsource d_fissionsource;

  //---------------------------------------------------------------------------
  // Smart pointers to other objects.

  /// State.
  SP_state d_state;

  /// Mesh.
  SP_mesh d_mesh;

  /// Angular mesh.
  SP_angularmesh d_angularmesh;

  /// Materials.
  SP_materials d_materials;

  /// M operator.
  SP_moment_to_discrete d_m_operator;

  /// Moments indexer.
  SP_momentsindex d_momentsindex;

  //---------------------------------------------------------------------------
  // >>> Implementation

};

} // end namespace slabtran

//---------------------------------------------------------------------------//
// INLINE FUNCTIONS
//---------------------------------------------------------------------------//

//#include "tran/SweepSource.i.hh"

#endif /* SWEEPSOURCE_HH_ */

//---------------------------------------------------------------------------//
//              end of SweepSource.hh
//---------------------------------------------------------------------------//
